function neuronDensityCloud(mode)
% NEURONDENSITYCLOUD
%
% gs July 2004
% ******************************************

if nargin == 0
    mode = 'UV beam 1';
%     mode = 'UV beam 2';
%     mode = 'glutamate';
%     mode = 'neuron cloud';
%     mode = 'uncaging volume';
end

% figure
figure(227);
set(gcf, 'Position', [34   699   560   420], 'Color', 'w');
colormap(winter);

% axes
delete(get(gcf, 'Children'));
axes;
xVal = 400;
xlim = [-xVal xVal];
ylim = [-300 300];
zlim = [-500 200];
set(gca, 'XLim', xlim, 'YLim', ylim, 'ZLim', zlim);
hold on;

% top and bottom of slice
X = [-xVal -xVal xVal xVal];
Y = [300 -300 -300 300];
Z = [0 0 0 0];
C = 'y';
htop = fill3(X,Y,Z,C);
set(htop, 'FaceColor', 'y', 'FaceAlpha', 0.2);
%
X = [-xVal -xVal xVal xVal];
Y = [300 -300 -300 300];
Z = [-300 -300 -300 -300];
C = 'y';
hbott = fill3(X,Y,Z,C);
set(hbott, 'FaceColor', 'y', 'FaceAlpha', 0.3);

% laminar boundaries


switch mode
    case 'UV beam 1'
        % hourglass cylinder -- UV beam
        x = -2:.01:2;
        n = 100;
        [X,Y,Z] = cylinder(x.^2+2,n);
        Xoffset = 0;
        Yoffset = 0;
        Zoffset = 0;
        X = X * 2 + Xoffset;
        Y = Y * 2 + Yoffset;
        Z = Z - 0.5;    % center in Z
        Z = Z * 1000 - 0;      % expand and offset in Z; focus depth in um into the slice
        hcyl = surf(X,Y,Z);
        set(hcyl,'LineStyle', 'none', 'FaceLighting','gouraud','FaceColor','interp',...
            'FaceAlpha', 0.5, 'AmbientStrength',0.5);
        
    case 'UV beam 2'
        % conical -- UV beam
        x = -2:.01:2;
        n = 100;
        [X,Y,Z] = cylinder(x.^2+2,n);
        X = X * 2 + -400;
        Y = Y * 2 + 0;
        Z = Z - 0.5;    % center in Z
        Z = Z * 1000 - 0;      % expand and offset in Z; focus depth in um into the slice
        hcyl = surf(X,Y,Z);
        set(hcyl,'LineStyle', 'none', 'FaceLighting','gouraud','FaceColor','interp',...
            'FaceAlpha', 0.5, 'AmbientStrength',0.5);
        
    case 'glutamate'
        % ellipsoid -- glutamate cloud
        xc = 0; % center position
        yc = 0;
        zc = -75;
        xr = 10; % radii
        yr = 10;
        zr = 75;
        n = 50;
        [x,y,z] = ellipsoid(xc,yc,zc,xr,yr,zr,n);
        Xoffset = 0;
        Yoffset = 0;
        Zoffset = 0;
        x = x - Xoffset;
        hellglu = surf(x,y,z);
        set(hellglu,'LineStyle', 'none', 'FaceLighting','gouraud','FaceColor','interp',...
            'FaceAlpha', 0.5, 'AmbientStrength',0.5);
        
    case 'neuron cloud'
        % large cloud of neurons
        sliceThickness = 300;
        numNeur = 500;
        rand('state',sum(100*clock));
        r = rand(numNeur, 6);
        for n = 1:numNeur
            xc = xlim(1) + (xlim(2)-xlim(1)) * r(n,1); % left position
            yc = ylim(1) + (ylim(2)-ylim(1)) * r(n,2);
            zc = 0 - sliceThickness * r(n,3);
            xr = 10; % radii
            yr = 10;
            zr = 10;
            n = 50;
            [x,y,z] = ellipsoid(xc,yc,zc,xr,yr,zr,n);
            hsomata = surf(x,y,z);
            set(hsomata,'LineStyle', 'none', 'FaceLighting','gouraud','FaceColor','interp',...
                'AmbientStrength',0.5);
            if abs(xc) < 25 & abs(yc) < 25 & abs(zc) < 25
                set(hsomata, 'FaceAlpha', 1);
            elseif abs(xc) < 50 & abs(yc) < 50 & abs(zc) < 50
                set(hsomata, 'FaceAlpha', 0.8);
            elseif abs(xc) < 75 & abs(yc) < 75 & abs(zc) < 75
                set(hsomata, 'FaceAlpha', 0.6);
            elseif abs(xc) < 100 & abs(yc) < 100 & abs(zc) < 100
                set(hsomata, 'FaceAlpha', 0.4);
            elseif abs(xc) < 150 & abs(yc) < 150 & abs(zc) < 150
                set(hsomata, 'FaceAlpha', 0.2);
            else
                set(hsomata, 'FaceAlpha', 0.05);
            end
%             if abs(xc) < 150 & abs(yc) < 150 & abs(zc) < 150
%                 if abs(xc) < 100 & abs(yc) < 100 & abs(zc) < 100
%                     if abs(xc) < 75 & abs(yc) < 75 & abs(zc) < 75
%                         if abs(xc) < 50 & abs(yc) < 50 & abs(zc) < 50
%                             if abs(xc) < 25 & abs(yc) < 25 & abs(zc) < 25
%                                 set(hsomata, 'FaceAlpha', 1);
%                             else
%                                 set(hsomata, 'FaceAlpha', 0.8);
%                             end
%                         else
%                             set(hsomata, 'FaceAlpha', 0.6);
%                         end
%                     else
%                         set(hsomata, 'FaceAlpha', 0.4);
%                     end
%                 else
%                     set(hsomata, 'FaceAlpha', 0.2);
%                 end
%             else
%                 set(hsomata, 'FaceAlpha', 0.05);
%             end
        end
        
    case 'uncaging volume'
        % ellipsoid -- uncaging volume
        xc = 0; % center position
        yc = 0;
        zc = -75;
        xr = 50; % radii
        yr = 50;
        zr = 75;
        n = 50;
        [x,y,z] = ellipsoid(xc,yc,zc,xr,yr,zr,n);
        Xoffset = 0;
        Yoffset = 0;
        Zoffset = 0;
        x = x + Xoffset;
        hell = surf(x,y,z);
        set(hell,'LineStyle', 'none', 'FaceLighting','gouraud','FaceColor','interp',...
            'FaceAlpha', 0.5, 'AmbientStrength',0.5);
end

% lighting
light('Position',[0 -1 1],'Style','infinite');

% axes
set(gca, 'XLim', [-xVal xVal], 'YLim', [-300 300], 'ZLim', [-500 200], ...
    'CLim', [-125 25], 'Projection', 'perspective');
set(gca, 'CameraPosition', [-583.438 -3147.95 -150], 'CameraViewAngle', [8.76548], 'View', [-10.5 0]);
daspect([1 1 1]);

box on ;

xlabel('vertical distance');
ylabel('horizontal distance');
zlabel('depth in slice');

